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ABSTRACT 

We revisit the linear MRI in a cylindrical model of an accretion disk and uncover 
a number of attractive results overlooked in previous treatments. In particular, we 
elucidate the connection between local axisymmetric modes and global modes, and 
show that a local channel flow corresponds to the evanescent part of a global mode. 
In addition, we find that the global problem reproduces the local dispersion relation 
without approximation, a result that helps explain the success the local analysis enjoys 
in predicting global growth rates. MRI channel flows are nonlinear solutions to the 
governing equations in the local shearing box. However, only a small subset of MRI 
modes share the same property in global disk models, providing further evidence that 
the prominence of channels in local boxes is artificial. Finally, we verify our results via 
direct numerical simulations with the Godunov code RAMSES. 
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1 INTRODUCTION 

The magnetorotational instability (MRI) remains the prin¬ 
cipal mechanism facilitating turbulence, and consequently 
mass accretion, in astrophysical disks. Since the seminal pa¬ 
per of Balbus and Hawley (1991, hereafter BH91), numerical 
simulations of MRI turbulence have become increasingly so¬ 
phisticated and comprehensive; now computations in global 
disk geometry are relatively common (e.g. Penna et al. 2010, 
Hawley et al. 2013, Parkin Sz Bicknell 2013) as is the inclu¬ 
sion of a panoply of physical effects (e.g. Jiang et al. 2013, 
Hirose et al. 2014, Lesur et al. 2014, Bai 2014). Despite this 
increase in complexity, BH91’s local and incompressible lin¬ 
ear theory continues to guide the interpretation of the sim¬ 
ulation results with surprising success (Flock et al. 2010, 
Hawley et al. 2011, Okuzumi & Hirose 2011). 

When it first appeared, however, the analysis of BH91 
provoked a debate concerning the local approach’s valid¬ 
ity when describing phenomena in the global geometry of 
an accretion disk. The emphasis of the debate focussed, 
in particular, on the small imaginary part exhibited by 
the MRI growth rates (neglected in the local analysis) as 
well as the boundary conditions’ impact on the magni¬ 
tude of the growth (also neglected; Knobloch 1992, Dubrulle 
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& Knobloch 1993). Ultimately, these concerns were shown 
to be somewhat exaggerated. Analyses in cylindrical mod¬ 
els, vertically stratified boxes, and general geometries indi¬ 
cated that the stability criterion was unaltered and the local 
growth rates good approximations in most instances (Ku¬ 
mar et al. 1994, Papaloizou & Szuszkiewicz 1994, Gammie 
& Balbus 1994, Curry et al. 1994, Curry & Pudritz 1995, 
Terquem Sz Papaloizou 1996, Ogilvie 1998). 

One issue that remains undeveloped is the question of 
how the BH91 local modes actually manifest in global disks. 
In other words: how do the local and global formalisms join 
up? This is especially important for the MRI channel modes, 
the fastest growing, because in the local analysis they exhibit 
no radial variation at all. This implies that their radial vari¬ 
ation is in fact global. What might their radial structure be 
— and does it matter? A second issue, also unexplored, is 
the nature of the modes that participate in simulations of 
global MHD turbulence. Simulated modes are never strictly 
local because of resolution constraints. We may ask how well 
such modes can be understood via the local approach. We 
may also want to construct a taxonomy of global modes and 
determine how they control the turbulence and other global 
aspects of the disk evolution (zonal flows, winds, magnetic 
flux diffusion, etc). 

In this paper we make explicit the connection between 
the local and global MRI, by revisiting the incompressible 
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axisymmetric cylindrical disk models employed by Dubrulle 
V Knobloch (1993), Coleman et al. (1994) and Curry et 
al. (1994). In so doing we uncover a number of attractive and 
overlooked results. First, we show that the global problem 
reproduces the (discretised) local dispersion relation with¬ 
out approximation. The boundary conditions only determine 
the discretisation, they do not generally influence the shape 
of the dispersion relation. Second, global MRI modes ex¬ 
tend over a limited range of disk radii before decaying in 
an outer evanescent region in which the differential rota¬ 
tion is too weak. We make clear that local channel modes 
can be identified with these evanescent portions, while local 
radially varying modes can be identified with parts of the 
same modes at smaller radii. Third, convenient analytic ap¬ 
proximations to the global growth rates are available via a 
matched WKBJ procedure. Fourth, in general, global MRI 
modes do not possess the nonlinear property of familiar local 
channel flows (Goodman V Xu 1994). Only modes localised 
to the inner edge of the disk remain acceptable solutions 
to the full equations when possessing nonlinear amplitudes. 
These results are verified numerically with a small set of sim¬ 
ulations using the Godunov code, RAMSES (Teyssier 2002, 
Fromang et al. 2006). Finally, we discuss the issues and ques¬ 
tions they raise, and future work that could begin to address 
them. 


2 LINEAR STABILITY ANALYSIS 

Within this paper we consider an annular slice through the 
midplane of the disk. The slice’s vertical thickness is taken to 
be much less than the disk’s vertical scale height. The back¬ 
ground vertical structure is hence neglected, and the model 
is quasi-global: ‘local’ in the vertical and ‘global’ in the ra¬ 
dial. For convenience we also dispense with the disk’s radial 
structure and treat the ionised fluid as incompressible. The 
resulting formalism remains mathematically tractable while 
retaining important global effects (boundary conditions and 
curvature). 

The first linear calculations in this set-up were under¬ 
taken by Velikhov (1959) and Chandrasekhar (1961), but it 
was not until the 1990s that accretion disks were explicitly 
modelled in this way (Knobloch 1992, Dubrulle V Knobloch 
1993, Kumar et al. 1994, Curry et al. 1994), the latter work 
generally confirming the local theory of BH91. More re¬ 
cently, Kersale et al. (2004) generalised the set-up to permit 
a radial flow of material through the inner boundary; and, 
while the classical MRI is recovered, spurious ‘wall modes’ 
are generated by the inner boundary condition. These modes 
control to some extent the ensuing nonlinear dynamics (Ker¬ 
sale et al. 2006). Finally, we note the work of Rosin V Mes- 
tel (2012) which included the Braginskii stress, and hence 
could describe the onset of instability in the weakly colli- 
sional plasma of the Galactic disk. 

In this section we first exhibit the main equations and 
repeat the classical local axisymmetric stability analysis for 
reference. Next the linearised equations in cylindrical geom¬ 
etry are reduced to a single second-order Sturm-Liouville 
equation, through which we display the similarities with the 


local problem. Numerical and analytic solutions are then 
derived. 

2.1 Governing equations 

We work with the equations of ideal incompressible MHD, 

a t u + u -Vu =-V$--VP t +d-B -VB, (1) 
p 47 Tp 

ftB + u • VB = B • Vu, (2) 

V • u = 0 (3) 

V • B = 0. (4) 

Here velocity and magnetic field are denoted by u and B 
respectively, p is the constant density, 4> is the gravitational 
potential of the central object, and Pt is the combined gas 
and magnetic pressure. The accretion disk is a circular annu¬ 
lus, with gas inhabiting cylindrical radii between r = ro and 
r = ri with ro < q. The vertical extent of the domain we 
set to positive and negative infinity. This setup describes mo¬ 
tions that possess short vertical lengthscales (shorter than 
the scale height) and up to long radial lengthscales (> ro). 

2.1.1 Equilibrium 

The governing equations admit the steady equilibrium solu¬ 
tion: 

u = rQ(r) e<£, B = B 0 e z , P t = P t (r), (5) 
The rotation frequency is a power law, 

O = O 0 (“) 

If q = 3/2 the disk is Keplerian and the background pressure 
is constant (no pressure gradient is required for radial force 
balance). When q ^ 3/2 the pressure gradient is non-zero, 
but its details are not important in what follows (see Curry 
et al. 1994 for further specifics). 

2.1.2 Perturbations 

The equilibrium is disturbed by axisymmetric modes of the 
type oc F(r)e lkzZ+st , where F is an r-dependent perturba¬ 
tion amplitude, k z is the (real) vertical wavenumber, and s 
is the (potentially complex) growth rate. The assumption of 


vertical locality means that k z ro 1. 

The ensuing linearised equations are 

su r — 2 flu# = —d r h' + ik z v 2 A b' r , (6) 

su,$ + (2 - q)Qu' r = ikzVAb'j,, (7) 

su z = —ik z b! + ik z v\b' z , (8) 

sb' r = ik z u r , (9) 

sb# = ik z u ,^ — qQb' r , (10) 

sb' z = ik z u z , (11) 

(l/r)d r (ru r ) + ik z u z = 0. (12) 

Here a prime indicates the perturbation, while b ; = RoB' 


and is hence dimensionless. The Alfven speed is defined by 
v\ — Bq/(47tp) and the enthalpy by h' = Pt/p • Note that 



these perturbations automatically satisfy the solenoidal con¬ 
dition on bk 

To complete the problem we must supply two boundary 
conditions, applied to the radial velocity u' r . The boundaries 
can be treated as hard walls or as stress free (Dubrulle & 
Knobloch 1993), in which case either u r or d r u' r is zero at 
r = 7*0 and r = n, or they may be treated as free surfaces 
(Curry et al. 1994), in which case a linear combination of 
u' r and d r u r is zero. As we will see later, it is not terribly 
important which we choose, only that the conditions are 
homogeneous. 

Note that the Alfven frequency k z VA associated with 
each mode is constant throughout the disk. In contrast, the 
orbital frequency Q(r) decreases with radius. Consequently, 
at sufficiently large radius magnetic tension dominates and 
the conditions for MRI become unfavourable. We expect a 
growing MRI mode of given k z to avoid such radii and in¬ 
stead emerge closer to the central body. If, however, the 
background magnetic field decays with radius faster than H, 
this need not be the case. 


2.2 Local axisymmetric dispersion relation 

In order to examine local modes, we choose a point r* and 
examine the behaviour of the gas in its immediate vicin¬ 
ity. For axisymmetric disturbances the orbital frequency Q 
can be regarded as constant in Eqs (6)-(12), and if their 
radial variation is small-scale then the cylindrical term in 
Eq. (12) may be dropped. This permits us to decompose 
the disturbances in Fourier modes oc e lkxr , and the local 
MRI dispersion relation for such modes is straightforward 
to derive: 

s 4 + [2 k 2 v\ + 2(2 — q)c 2 Ot\s 2 + k 2 v 2 A (k 2 v 2 A — 2ge 2 fl*) = 0. 

(13) 

Here fl* = n(r*), the rotation rate at the radius in which 
we’re interested, and 

£= (i+W (14) 

In Eq. (13) the radial wavenumber k x appears solely in 
e and then only in the ratio k x /k z . Note also that instances 
of e occur exclusively as factors of fl*. The familiar channel 
flows are obtained when k x /k z —> 0, with the modes exhibit¬ 
ing little (to no) relative radial variation. In this case, e —> 1 
and the fastest growth rate is achieved, qQ*/ 2. Because eQ* 
sets the timescale of the growth rate, modes with non-zero 
k x grow slower, as they possess e < 1. We call radially vary¬ 
ing modes ‘radial modes’ to distinguish them from channel 
modes. 

The local dispersion relation (13) exhibits an interesting 
connection between radial variation , on one hand, and radial 
location , on the other. By varying k x (and keeping k z fixed) 
we can ‘rescale’ the MRI timescale efl*, because the factor e 
depends on k x . A striking consequence is that a radial mode 
of a given e located at one radius n, associated with an 
orbital frequency of Q(r*) = fl*, possesses the same growth 
rate as a channel mode located at a different larger radius 
r Q > n, associated with an orbital frequency of Q(r 0 ) = eQ*. 



Radially-varying mode 

Figure 1 . Simple illustrations of a channel mode (upper panel) 
and a radial mode (lower panel). Blue circles represent fluid blobs, 
the black arrows indicate the initially vertical magnetic field. Red 
arrows show the direction of the fluid motion. 


As a consequence, we can ‘map’ modes of different k x , but 
different radial locations, onto one another. 

Physically, one can understand the radial modes by re¬ 
turning to the classical cartoon of the MRI mechanism (for 
e.g. BH91). The upper panel of Fig. 1 illustrates a channel 
mode, with the blue circles representing fluid blobs and the 
black lines indicating the magnetic field (initially vertical). 
Vertically varying radial perturbations lead to the periodic 
stretching of the magnetic field and, consequently, to angu¬ 
lar momentum exchange between the tethered blobs. The 
instability proceeds from the counterstreaming flows that 
ensue. 

However, if the mode possesses radial structure there 
will exist nodes in radius where the streams converge and 
diverge, illustrated in the lower panel of Fig. 1. At these 
radii the mode exhibits vertical deflection and pressure per¬ 
turbations which are not present in pure channel flow. Need¬ 
less to say, the development of vertical circulations impedes 
the instability mechanism exemplified in pure channel flow. 
During their vertical excursions fluid blobs stop extracting 
free energy because they are no longer exchanging angular 
momentum effectively. 


2.3 Global eigenvalue problem 


We now return to the linearised global equations (6)-(12) 
which can be reworked into a single second-order equation 
for the variable U = r 1//2 iv- Rescaling space by ro we find 


d 2 U 

dr 2 




U, 


(15) 



which should be solved on the domain r G [1, ri/ro\. Here 
the various parameters appearing in the problem, in addition 
to the unknown growth rate s, have been packaged into the 
convenient quantity £ defined via 


£ 


-2 


= 20q 


(q - 2)s 2 + q/tjvA 
(s 2 + klv 2 A ) 2 


(16) 


In the above, flo = U(ro), the rotation rate at the inner 
boundary. 

With homogeneous boundary conditions (such as im¬ 
penetrable walls or a free surface), Eq. (15) is in Sturm- 
Liouville form. The weight function is —k 2 r~ 2q and the 
eigenvalue is £~ 2 . Because the problem is Sturm-Liouville, 
we are assured of a discrete set of real eigenvalues e~ 2 which 
we order so that £o > £i > £2 > .... Once these are com¬ 
puted the associated growth rates can be obtained. Hence 
the problem is broken down into two steps: first calculate 
U and £ from (15), then calculate the growth rate s from 
(16). Note that the eigenfunction structure U depends only 
on k z ; it is oblivious to the magnetic field strength va unless 
it appears in the boundary conditions. 


2.3.1 Global dispersion relation 

Equation (16), which defines the eigenvalue e, can be re¬ 
worked into 


2.3.2 Stability criterion 

Next we rearrange Eq. (15) so that it resembles a 
Schrodinger equation, 

^+klf{r,e)U = 0, (18) 

where 

= " 4 r 2 fc 2 “ L ( 19 ) 

The problem now resembles that of a particle in a potential 
well, with the potential proportional to /. It is easy to show 
that / has an extremum at 

r=(f qk 2 /e 2 f mq - 1)] . (20) 

In order for there to be trapped waves, i.e. normal modes, 
the extremum must lie inside the disk. For simplicity we 
let its outer boundary go to infinity, then this condition, 
combined with (17), gives us a general instability criterion 
in terms of the vertical field 

> ;\/!' <21) 

where Ma is the Alfvenic Mach number of the background 
state, defined to be equal to Qoro/vA- The criterion (21) 
gives an upper bound on the magnetic field threading the 
disk. Because in most contexts the field is assumed to be 
subthermal and the disk assumed thin, we have Ma 1 and 
the criterion is automatically satisfied. A far more restrictive 
condition on the magnetic field issues from the disk’s vertical 
thickness (see BH91 and Gammie & Balbus 1994). 


s 4 + [2 k 2 v\ + 2(2 — q)£ 2 a Q^\s 2 + k 2 v\ (k 2 v\ — 2 ^£^Qq) — 0 , 

(17) 

which is almost identical to Eq. (13)! Remarkably, the ax- 
isymmetric global problem yields a variant of the local ax- 
isymmetric dispersion relation. The basic structure of the lo¬ 
cal problem is retained independently of the global specifics: 
curvature and boundary conditions. 

There are two differences. Instead of the function e, 
which depends smoothly on the continuous radial wavenum¬ 
ber k Xl the global problem possesses the discrete function 
s n which depends on the radial quantum number n. Un¬ 
like the local problem, we do not know how e n depends 
on n a priori; this information must be extracted from 
the differential equation (15) and depends on k z , q , and 
the boundary conditions. However, it is easy to show that 
e n < 1 (see Appendix A). Hence we order the eigenvalues 
asl > £0 > £1 > £2 > ••• >0. The second difference is 
that Q is set to Uo, the orbital frequency at the inner radius 
of the disk. This fixes the minimum timescale over which 
appreciable growth can happen. 

It should be emphasised that, in contrast to the claims 
of previous authors, the only role the boundary conditions 
play is to help set the discretisation of the dispersion relation 
(17). It does not alter the shape of the relation itself. That 
fundamental shape is determined by the local physics, as 
first explained in BH91. Moreover, as we show later, the 
larger rok z the more closely spaced the £ n and hence the 
less important the discretisation. 


2.3.3 Turning points and global and local modes 

At sufficiently large r the function / will be dominated by 
the last term in Eq. (19). In this region U ~ e~ kzV and the 
mode is evanescent. Physically this makes sense: sufficiently 
far out in the disk the differential rotation is too weak to 
compete with magnetic tension and the MRI mechanism 
is suppressed. Active MRI modes will shun such a region 
and instead localise at smaller radii where the conditions 
for instability are more favourable. In Fig. 2 we present a 
representative eigenfunction illustrating this morphology. 

The boundary between the evanescent region and the 
‘MRI region’ is given by the turning point of Eq. (18), i.e. 
the radius at which f — 0, denoted by r tp . For general q 
and k z this must be determined numerically. However, in 
the plausible limit of k z large, 

rtp » £~ 1,q . (22) 

In fact, this radius corresponds to an orbital frequency of 
Q = eflo- 

Combining this last piece of information with the local 
and global dispersion relations reveals an interesting corre¬ 
spondence. A global mode, associated with a given e n , shares 
the same dispersion relation (17) as that of a local channel 
mode situated at the global mode’s turning point (13). This 
is because the local orbital frequency Q* at this point is pre¬ 
cisely £ n Tlo. Consequently, we are invited to identify a local 
channel mode as the small section of a much larger global 




Figure 2. A higher order global eigenfunction of (15). The turn¬ 
ing point rt p is indicated as are regions of the global mode that 
we identify with local channel and radial modes. Note also the 
monotonic decrease in the local radial wavenumber k r with ra¬ 
dius. 

mode — the section near evanescence. Indeed, we may define 
a (radius dependent) radial wavenumber via k r « k z f 1 ^ 2 , 
and at the turning point k r = / 1//2 = 0, by definition, which 
is in accord with the channel mode’s lack of any radial struc¬ 
ture. 

Meanwhile, what of local radial modes? First, consider 
a radial mode located at r « r*o, i.e. very near the inner 
boundary, and possessing a k x so that e « e n , for the same 
n as above. This particular radial mode now shares the same 
dispersion relation as the previous channel mode, and hence 
the same dispersion relation as the global mode, Eq. (17). 
Hence we may also identify this local radial mode as part of 
the same parent global mode — but the part of that mode 
near the inner boundary. Moreover, we can do exactly the 
same with other radial modes of smaller k x . But for each 
different k x we must change the local mode’s radial location 
so that the timescale eQ is kept constant. 

This explains the connection between radial variation 
and radial location that appeared in the local dispersion re¬ 
lation of Section 2.2. This connection arises because a global 
mode of given n is constituted from a set of local modes of 
differing k Xl each at a different radial location. In Fig. 2 this 
idea is sketched out. 

It is worth stressing that the growth rate of the global 
mode is limited by the local physics at the mode’s periphery, 
at r — 7* t p. Though the mode extends over regions where the 
local growth rate can be faster, the mode can only grow as 
fast as its outer ‘edge’. Put another way: the structure, as a 
whole, can only grow as fast as its slowest component. This 
helps remove some of the ambiguity when attributing local 
growth rates to global simulations. It is clear that at any 
given radius, growth is controlled by the mode whose turning 
point falls at that radius. No other mode that extends to this 
location can grow faster (faster growers are localised to radii 
closer in). Moreover, as we see from Figs 2 and 3, modes’ 
amplitudes are maximal near r — r tp . Indeed, numerical 
simulations verify this. At any given radius, growth occurs 
at the rate of the local channel mode (Hawley 2001), which 


of course is also the rate of the global mode whose turning 
point r tp falls there. 


2.4 Global solutions 

2.4.1 Numerical solutions 

Having sketched out the background details, we present in 
this subsection a few numerical solutions to (15) subject to 
the hard wall boundary conditions, U — 0 at r = 1, n/r*o. 
Our domain is set to be r £ [1,10]. The only parameter that 
appears explicitly in (15) is the vertical wavenumber k z . We 
let it equal 10 for our main results. The numerical technique 
we use is a pseudo-spectral Chebyshev method (Boyd 2002), 
which approximates the differential equation by a matrix. Its 
eigenvalues may be obtained by the QZ algorithm (Golub & 
van Loan 1996). The boundary conditions are encoded in 
the matrix via boundary bordering. 

In Fig. 3 appear the first ten eigenfunctions for k z = 10. 
The associated e n for the first four are 

{0.6413, 0.4976, 0.4169, 0.3628}, 
with corresponding turning points r*t p 

{1.345, 1.593, 1.792, 1.966}, 

where the mode transitions to an evanescent wave. As antic¬ 
ipated, larger n correspond to modes that extend over more 
of the domain. 

In order to compute explicit growth rates we need to 
specify the strength of the vertical magnetic field, this can 
be measured in terms of either the dimensionless combina¬ 
tion VA{k z /Q) or the Alfvenic Mach number Ma • Selecting 
the first option, we plot in Fig. 4 the growth rates of the 
first four modes as functions of va and fixed k z = 10. These 
are just four copies of the standard local MRI relation. The 
different curves correspond, not to differences in the modes’ 
vertical structure — this is held fixed — but to different 
radial structures. The greater the radial quantum number n 
the more radial structure, but most importantly the greater 
r tp and hence the slower the dynamical timescale. This ex¬ 
plains why higher n modes grow slower. It should be appre¬ 
ciated that the fastest growing mode need not determine the 
evolution of the disk globally. The fastest growing modes are 
localised at small radius. Larger radii will be driven by the 
slower modes that extend to them. 

Lastly, we fix va and see how the growth rates vary as a 
function of k z . In Fig. 5 we plot s as a function of k z , different 
n, and for a weak field fixed by Ma = 50. These curves are 
similar but not exactly those of the standard MRI, seen in 
Fig. 4. That is, we cannot simply rescale k z ro by k z VA/& 0 
and recover the same relation. This is because £ n — e n (k z ). 


2.4.2 Special solutions 

For certain special cases, Eq. (15) can be solved exactly. 
When the disk exhibits a rotation profile of q — 1, roughly 
similar to the Galaxy, we find that 

U = r 1/2 K„{k z r), with v = ^1 - fcf/e 2 , (23) 






































Figure 3. The first ten MRI eigenfunctions when k z = 10. In each case, the modes are oscillatory within the region 1 < r < rt p and 
clearly evanescent outside these. Real parts are plotted; the imaginary parts are zero. 





Figure 4. MRI growth rates for five different radial quantum 
numbers and when k z = 10 but va, the strength of the back¬ 
ground magnetic field, varies. Note that va has been scaled using 
a fixed k z ; this is so the classical MRI dispersion curves are easier 
to see. 


and where K u (x) is the modified Bessel function of the sec¬ 
ond kind. In obtaining (23) we have taken the outer bound¬ 
ary to infinity and then applied the outer boundary condi¬ 
tion, u r — r~ 1 / 2 U —>> 0. The eigenvalue equation for e is 
K„(k z ) = 0 in the case of a hard inner wall. These solu¬ 
tions were explored first by Dubrulle & Knobloch (1993), 
and more recently by Rosin &; Mestel (2012). 

In the limit of large k z , approximations to the eigenval¬ 



Figure 5. MRI growth rates for the first 10 radial quantum 
number n when va is fixed and k z is varied. We have set 

M a = ro^o/vA = 50. 

ues can be obtained from 

£n«l + a n 2“ 1/3 fcj 2/3 , (24) 

where a n is the n’th negative root of the Airy function Ai(a;). 
A brief derivation of this expression is given in Appendix A2 
(see also Rosin & Mestel 2012). Note that in this limit the 
spacing between neighbouring eigenvalues becomes tiny, and 
thus the dispersion relation (17) is effectively continuous: 
local physics dominates the global mode. 









































































When the rotation profile is Keplerian, q — 3/2, there 
is no general analytic solution to Eq. (15). However, in the 
limit of small k z we have 

U « r 1/2 J 2 (?L r “ 1/2 ) - (25) 

where J^{x) is the Bessel function of the first kind and sec¬ 
ond order. The ensuing eigenvalues are 

£n 2k z /b n , (26) 

where b n is the n’th root of J 2 (x). Of course, the limit of 
small k z is not the most appropriate for our cylindrical disk 
model, as it indicates the mode is global in the ^-direction. 
The result nevertheless offers a useful test of the numerical 
solver in Section 2.4.1. 


2-4-3 WKBJ solutions 


In this section we present asymptotic solutions to (15) in the 
limit of large vertical wavenumber k z and for general q. As 
explained earlier, this is a natural limit because we expect 
k z ro to be large. 

The solution is oscillatory between r — 1 and r — 
£~ x ^ q — rtp, and it is evanescent for r > n p . Within the 
former region we employ the standard WKBJ ansatz: 


u = I/I 1/4 cos 


k z [ \f\ 1/2 dr + ir/4 

J n P 


(27) 


The phase shift of 7r/4 comes from matching across the turn¬ 
ing point (see for e.g. Riley et al. 2006). To obtain the eigen¬ 
value equation for £ we next impose the hard-wall boundary 
condition at r — 1, i.e. U — 0. This leads to 



where n is the radial quantum number and in which we 
have set f(r) = ( er q )~ 2 — 1 to leading order in k z . Notice, 
that for large k z neighbouring n modes posses eigenvalue 
equations that differ only marginally on the right hand side. 
As a consequence, the values of the corresponding e n are 
extremely close to one another, as in Eq. (24). 

The integral in (28) can be written in terms of special 
functions by introducing the integration variable £ = s 2 r 2q . 
The eigenvalue equation then becomes 


2n q 


*(£-*,§)-*.■ (^-i§) = T^ 1 A >+i) 


(29) 

where B(x,y) and B^x, y ) are the complete and incomplete 
beta functions (Abramowitz V Stegun 1964). Though it is 
relatively straightforward to numerically solve this nonlinear 
equation, we have the following convenient approximations. 
The fastest growing modes possess 5 near 1, which may be 
approximated by 


e n ~ 1 - | [37 iq (n + |)] 2/3 k z 2/3 . (30) 

Note the similarity to Eq. (24 ) 1 . For small and intermediate 


1 Incidentally, when q = 1 this equivalence provides a neat an¬ 
alytic approximation to the zeros of Ai(x). For a more direct 
derivation see Fabijonas and Olver (1999). 


£ a different expansion of the beta functions yields 



1 7r 

2(1 + q) £ + k z 



= 0 , 


(31) 


where A is the number 



(32) 


and T[x] is the gamma function. A Keplerian rotation law 
reduces Eq. (31) to a quintic polynomial equation. 

For comparison with the computations in Section 2.4.1, 
the WKBJ eigenvalues £ are 


{0.645, 0.499, 0.418, 0.363}, 


which agree reasonably well. The agreement improves as k z 
increases. Values of e gathered from either (30) or (31) may 
be input into the dispersion relation (17) and the growth 
rates computed without the need to numerically solve the 
ODE, which is the main benefit of the WKBJ approach. 


3 ARE GLOBAL MODES NONLINEAR 
SOLUTIONS? 

In classical and vertically stratified shearing boxes, chan¬ 
nel flows are nonlinear solutions in the incompressible and 
anelastic regimes, respectively (Goodman V Xu 1994, Lat¬ 
ter et al. 2010). It is then natural to ask if global modes in 
cylindrical geometry possess an analogous property. Indeed, 
the initial stages of some nonlinear simulations (e.g. Haw¬ 
ley 2001) exhibit strong counterstreaming flows in the inner 
parts of the disk that suggest this might be the case. In this 
section, however, we show that only a small subset of the 
linear global modes can be said to be approximate nonlin¬ 
ear solutions, and furthermore these are localised very close 
to the inner boundary. 

The key feature of local channel flow is the strong sep¬ 
aration between their radial and vertical lengthscales. We 
hence introduce the small parameter 

S = l/(k z Xr), (33) 

where A r is the modes’ characteristic radial lengthscale. 
From Eqs (6)-(12) the following scalings can be derived 

/ c / / / 

u z ~ o u r , rsj u r 

b' z ~5b' r , b't-K, 

in addition to h! ~ S 2 (X r Q)u r . Thus in the limit of small 
5 the morphology of global modes resembles local channel 
flows: vertical velocity and magnetic field is minimised, as 
is the pressure perturbation. Perfect channel flows cannot 
be obtained, however, because of the cylindrical terms in 
conjunction with incompressibility and the solenoidal con¬ 
dition. Some (small) radial variation in the mode structure 
must give rise to (small) vertical motion or held. In the local 
approximation k z and A r can be chosen independently and 
5 can be made as small as required. This is not the case for 
global modes, however, and generally A r = A r (k z ). We de¬ 
fer estimates of the magnitude of 5 to later in this section; 



for the moment we assume that 5 can be made sufficiently 
small. 

Let us now inspect the size of the non-linearities as¬ 
sociated with the modal solutions. We want to estimate at 
what point the nonlinearities become important. Assuming 
the above scalings, the ratio of the nonlinear advective term 
in the radial component of the momentum equation to the 
linear terms is 


(u' • Vu') r 
su' r 


u r d r u r + u z ik z u r — u£/r 
Qu' r 





(34) 

(35) 


To progress further, we assume that the mode began growing 
at a fraction a of the background Alfven speed va, and thus 
|u'| oc ae st VA (see Goodman & Xu 1994). Placing this in 
the above scaling gives 


(iT • Vu / ) r 

su' r 


VA st 

x r n ae 

5 a e st , 


(36) 

(37) 


an estimate of e n in the limit of large k z , and we get A r ~ 
k z 2/ \ This returns 

(42) 


a rather weak scaling. As a consequence, small values of 5 
require exceptionally large values of k z . Indeed, only modes 
with k z > 10 3 possess anything resembling the nonlinear 
property seen in local boxes. Moreover the radial extent of 
such modes scale like k z 2 ^ 3 and so they are effectively lo¬ 
calised to the inner edge of the disk. As a consequence, the 
vast majority of an accretion disk will never experience the 
nonlinear property of the linear MRI modes. 

We finish by discussing alternative mechanisms that dis¬ 
rupt those few global modes that do possess an approximate 
form of the nonlinear property. Compressibility is perhaps 
the primary mechanism in local boxes (Latter et al. 2010). 
We can estimate when compressibility becomes important 
by setting |u'| ~ HQ. This gives a disruption time of 





(43) 


where we have assumed that VAk z /Q < 1, in order for the 
MRI to work. A similar argument gives the same scaling for 
the other components of the advective term, as well as the 
B • VB and mixed terms. For instance, using (6)-(12), 


(B' • VB') r /(47rp) 2 KdrK + b' z \k z b' r - b$/r 

su' r Qu' r 

v\k 2 u r 

~ n 2 nv 

u r 


(38) 

(39) 

(40) 


and the same scaling as (37) is recovered. 

In all cases the quadratic nonlinearities go as Sae st rel¬ 
ative to the linear terms. This immediately gives us a con¬ 
dition for when they are important and when the MRI de¬ 
viates from the linear channel structure. We set Sae st to 1, 
and compute the time at when this occurs. The critical time 
£1 is 


Compressibility intervenes before the quadratic nonlineari¬ 
ties only when £2 < £ 1 , which occurs if H < A r , a possible 
regime for modes of moderate to small k z . 

Alternatively, a runaway channel may be destroyed 
through the action of a parasitic mode feeding off its strong 
shear and magnetic energy (Goodman & Xu 1994, Pessah 
& Goodman 2009, Latter et al. 2009, 2010). By analogy 
with local boxes, we expect a variant of the vertical Kelvin- 
Helmholtz instability to be the fastest growing parasite. 
However, as argued in Appendix B in Latter et al. (2010), 
the orbital shear dramatically weakens the ability of the 
non-axisymmetric parasites to successfully destroy a chan¬ 
nel mode: usually, a channel is safe to grow to equipartition 
strengths. In cylindrical geometry, it is less clear whether 
parasitic modes are equally ineffective. 


4 NUMERICAL SIMULATIONS 



It follows that the nonlinearities only become important 
once the mode amplitudes have grown a factor 5 _1 greater 
than the background field. Thus the smaller 5 the more dra¬ 
matic the amplification of the linear cylindrical modes, and 
the deeper they penetrate the nonlinear regime. 

In fact, for most parameter choices we find only 5 < 1. 
This is because the larger we take the value of k z , the smaller 
the corresponding A r , and as a result 5 need not be small. 
Most modes, consequently, cannot be said to possess the 
nonlinear property of MRI channel flows: once their ampli¬ 
tudes reach that of the background field, quadratic nonlin¬ 
earities become important and the modes break down. 

However, for small n and exceedingly large k z the radial 
and vertical scales can be separated and <5 1 is achievable. 

The scaling of S with large k z is straightforward to obtain. 
First assume that for small n the radial scale of the mode is 
roughly equal to the distance between r — 1 and its turning 
point r = 7* t p. Thus A r & — 1. Equation (30) provides 


We now present a series of numerical simulations that il¬ 
lustrate some of the properties discussed above. The first 
objective is to follow the growth of the normal modes cal¬ 
culated in Section 2, and then study the breakdown of the 
linear regime. We start by describing the numerical method 
and the setup of our simulations. 

4.1 Method and setup 

We use the finite volume code RAMSES (Teyssier 2002, Fro- 
mang et al. 2006) to solve the axisymmetric compressible 
MHD equations in a 2D cylindrical coordinate system (r, z). 
Compressibility remains small during the linear phase of the 
mode evolution, which RAMSES accurately captures at the 
cost of a smaller timestep. In addition, finite volume codes 
are now widely used in studying the properties and con¬ 
sequences of the MRI in astrophysical disks; our use of a 
similar code will thus ease connections with previously pub¬ 
lished results and aid future work. 

The strategy we apply is the following: take a disk model 
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Figure 6. Color contours showing the spatial distribution in the 
( r,z ) plane of B r /Bo for the mode with k z = 40 and n = 5 
at time t = 2. The contour lines overplot B r /Bo at time t = 0 
and show that the mode has grown unperturbed over that period. 
Contours are drawn from —6 X 10 -5 to +6 x 10 -5 every 2 x 10 -5 . 
Negative contours are dashed and the zero contour is omitted. 



in dynamical equilibrium, superpose the eigenmodes calcu¬ 
lated in Section 2.4.1 with a very small amplitude, follow 
their growth and compute the associated growth rate. We 
take the simplest possible initial disk state: an ideal gas 
of uniform density, uniform temperature (or, equivalently, 
speed of sound Co), and adiabatic index of 5/3, in Keplerian 
rotation around a central point mass M. Units are fixed by 
choosing GM — 1, while the location of the grid’s inner 
boundary is at r = ro = 1. We use Co — 0.1 though the 
remainder of this section. 

Even though the setup is fairly standard, two aspects 
deserve a more detailed discussion: the nature of the ther¬ 
modynamic perturbations and the boundary conditions. The 
eigenmodes described earlier are in the incompressible limit 
and involve a non-zero pressure perturbation (so that the 
velocity divergence associated with the eigenmode remains 
zero). However, in compressible simulations, associated den¬ 
sity perturbation will feed into the momentum equation and 
create additional, non modal velocities perturbations, com¬ 
plicating the interpretation of our results. To avoid that 
problem, we introduce temperature perturbations that per¬ 
mit non-zero pressure perturbations but which keep the den¬ 
sity constant and uniform. 

The boundary conditions we used are periodic in the 
vertical direction. This is possible in our setup because we 
neglect the density’s vertical stratification (see Section 2.1). 
The radial boundary conditions are more subtle to imple¬ 
ment. As discussed earlier, the eigenmodes are confined to 
the inner disk. Thus, at the outer radial boundary, we force 
the variables to take their equilibrium, unperturbed values. 
At the inner boundary, we decompose the variables in the 
ghost cells into the sum of their equilibrium values (known 
analytically) and a perturbation. The latter are chosen using 
the first active cells of the grid with the constraint that they 
should satify the symmetry of the eigenmodes: zero-gradient 
for the vertical velocity and magnetic field perturbations and 
anti-symmetric for the radial and azimuthal components of 
the velocity and magnetic field (of course, the density pertur¬ 


bation vanishes). With these boundary conditions, the disk 
equilibrium is conserved to within machine accuracy and, as 
we shall see below, we can follow the MRI eigenmodes to 
saturation. 


4.2 Linear growth 


In this section, we numerically compute the growth rate 
of the normal modes for a vertical magnetic field whose 
strength is such that /3, the ratio between thermal and mag¬ 
netic pressure, is equal to 200 (or, equivalently, Ma = 100 
given our value of the sound speed). For a given eigenmode 
(defined by the values of k z and n), we compute the mode’s 
spatial structure and theoretical growth rate as described 
in Section 2.4.1. At t — 0, we add the perturbation associ¬ 
ated with that mode to the equilibrium disk structure and 
start the simulation. Its amplitude is such that the maxi¬ 
mum perturbed radial magnetic field amounts to 10 -4 Ho. 
The computational domain extents radially out to n = 3 
and the vertical size of the domain is set equal to the verti¬ 
cal wavelength of the mode. The value of r± is here chosen 
large enough to ensure that all variables reach their equilib¬ 
rium value (i.e the mode amplitude goes to zero) well within 
the domain. In order to evaluate the amplitude of the mode 
during a simulation, we compute the time evolution of the 
volume averaged rms of the radial component of the mag¬ 
netic field B r : 


C 2 (B r ) = 


/ ff rB^drdz 
V ffrdrdz 


1/2 


(44) 


The instantaneous growth rate a n is then defined as the time 
derivative of log(£ 2 )- We denote by dn its value averaged in 
time over the first two orbits of the simulation (this short 
timescale ensures that the system remains within the linear 
phase). 

To illustrate our general results we consider the spe¬ 
cific mode k z = 40 and n = 5. Its theoretical growth rate is 
a = 0.384. After simulating the mode’s evolution for various 
grid sizes, we found that a resolution of (7V r , N z ) = (256,16) 
is required to capture correctly the growth rate. In this case, 
ga — 0.380, which corresponds to the theoretical growth rate 
to better than a percent. In addition, the mode structure is 
not modified during the linear growth. This is illustrated 
in Fig. 6, which shows a good correspondance between the 
mode spatial structure in the (r, z) plane at t = 0 (contour 
lines) and t — 2 (color contours), despite an amplification 
by about three orders of magnitude. We note that the rela¬ 
tive density fluctuations remain smaller than 10 -4 at t — 2, 
which explains the good agreement between the compress¬ 
ible simulation and the incompressible linear analysis. De¬ 
creasing the number of vertical cells to 8 gives ga — 0.339, 
and the vertical structure of the mode is incorrectly de¬ 
scribed. Likewise, when decreasing the radial resolution by 
a factor of two, so that N r = 128, the mode structure in the 
inner parts of the disk is modified, even though oa — 0.371, 
which is a reasonable estimate. The mode’s radial wave¬ 
length decreases inward and requires a fine resolution to be 
properly captured. 

Based on these results, we kept a fixed resolution of 
(N r , N z ) — (256,16) and systematically evaluated ga for an 
































Figure 7. Growth rates of the unstable modes in the ( k z ,n ) plane for the case (3 = 200 and co = 0.1 determined using the linear analysis 
described in Section 2 ( left panel) and using an ensemble of 240 numerical simulations ( left panel) that have a common spatial resolution 
(. N r ,N z ) = (256,16). 



Figure 8. Instantaneous growth rates normalised by the theo¬ 
retical growth rate as functions of B r /Bo, the maximum value of 
the perturbed radial field. In each of the three cases n = 0, but 
k z = 90 ( blue curve), k z = 900 ( red curve), and k z = 9000 ( green 
curve). The corresponding plasma betas are (3 = 200, 2 x 10 4 , 
and 2 x 10 6 . 


ensemble of 240 simulations in which we varied the vertical 
wavenumber k z between 20 and 160 and the mode number 
n between 0 and 12. The results are summarized in Fig. 7 
and show excellent agreement between the theoretical ex¬ 
pectations (left panel) and growth rates estimated from the 
numerical simulations (right panel). Discrepancies occur for 
the slowest growing modes because they are typically over¬ 
whelmed by faster growing modes seeded by truncation er¬ 
rors. Except for these cases, we conclude that a finite volume 
compressible code can reproduce the results of the incom¬ 
pressible linear analysis. 


4.3 Nonlinear Saturation 

We now determine what happens when the modes enter 
the nonlinear regime. The simulations remain axisymmet- 
ric, even though realistically the final saturated state will be 
non-axisymmetric. Our aim, however, is not to characterize 
the turbulence that follows the linear mode breakdown, a 
formidable task that is well beyond the scope of this paper. 
Rather, we ask the simpler questions: what is the ampli¬ 
tude of a mode when it stops growing exponentially with 
time? Do there exist linear modes that penetrate the non¬ 
linear regime while maintaining their structure, as predicted 
in Section 3? 

We performed three simulations for which [3 — 200, 2 x 
10 4 and 2 x 10 6 . We focused on n — 0 and select in each case 
a mode with growth rate close to the maximum rate, thus 
k z = 90, 900 and 9000, respectively. Note that only in the 
last case is 5 unambigously small, being ~ 0.05 according to 
Eq. (42). The vertical and radial extent of the computational 
domain is changed in each of the three simulations to best 
accomodate each mode, as described in the previous section. 

Figure 8 displays the instantaneous growth rates a n nor¬ 
malized by the theoretical growth rates a as a function of 
Br/Bo , where B r denotes here the maximum of the radial 
magnetic field fluctuation. Because B r /Bo increases mono- 
tonically it may be regarded as a proxy for time. The non¬ 
linear regime corresponds to B r /Bo > 1. The figure clearly 
shows that the larger k z , the greater the amplitude achieved 
by the mode before it breaks down. This is consistent with 
the prediction of Section 3, which states that quadratic non- 
linearities intervene only once a mode grows to d -1 the back¬ 
ground. For the case k z = 9000, the estimated maximum 
amplitude of the perturbed field is 20, in fair agreement 
with Fig. 8, which shows that a n /a reaches 0.9 at this point. 
This mode, in particular, has penetrated significantly into 
the nonlinear regime. This is not so convincingly the case 
for the lower k z modes, as expected. 

In Fig. 9 we compare the distribution of the radial mag¬ 
netic field at t = 0 (top row) and at a later time when B r 















has grown to amplitudes larger than the background verti¬ 
cal field (bottom row) for k z = 90 and k z — 9000. The figure 
shows that the low k z mode has become strongly disturbed 
by t = 3.2, whereas the field of its large k z counterpart 
retains a structure close to linear, despite possessing an am¬ 
plitude 10 times greater than the background. 

Finally, we check the role of compressibility, which 
should be important when H < A r . In code units, H — 0.1 
at r = ro, whereas A r ~ 0.1 and A r ~ 0.001 for the two 
modes k z = 90 and k z = 9000. Compressibility is certainly 
subdominant in the breakdown of the large k z mode, which 
we attribute to quadratic nonlinearities. In the breakdown 
of the low k z , on the other hand, compressibility and non¬ 
linearity are of equal importance. 


5 CONCLUSION 

We have presented a linear analysis of the MRI in cylindri¬ 
cal geometry in order to make clear the connection between 
the global and the local theories. In particular, we show 
that each local channel mode corresponds to the evanescent 
part of a global mode. Local radially varying modes, on the 
other hand, correspond to sections of the same global mode 
at smaller radii. Moreover, we show that only global axisym- 
metric modes of extremely large vertical wavenumber k z are 
approximate nonlinear solutions to the governing equations. 
As these modes are localised to the inner boundary, most 
of the disk never experiences the ‘nonlinear property’ of the 
MRI. Direct simulations with RAMSES verify the last point, 
which also provides a useful numerical check on global codes 
generally. Our results raise a number of questions and issues, 
which we now list. 

First, one of the most notable features of channel flows 
in local boxes is that they are nonlinear solutions to the 
governing equations. But as this property fails to appear in 
global disk models it is clear that this feature is an artefact 
of the local approximation. A natural question is then: how 
seriously does this artefact distort the nonlinear dynamics 
of the MRI in shearing boxes? It is well known that channel 
flows strongly influence the MRI saturation in certain cir¬ 
cumstances (e.g. Sano V Inutsuka. 2001, Bodo et al. 2008, 
Lesaffre et al. 2009, Murphy & Pessah 2015). How repre¬ 
sentative of global MRI turbulence is local MRI turbulence, 
given the artifical prominence of channel flows in the latter? 

Second, how non-local is the onset of the MRI and its 
ensuing turbulence in global disks? At a given radius r*, the 
fastest growing disturbance is the normal mode that has its 
turning point there. But such a mode extends from the in¬ 
ner boundary r = 1 to r = r*, so its growth, and subsequent 
behaviour, may be influenced by all the shorter time-scale 
activity on the intervening shorter radii. Indeed, because 
different global modes encompass overlapping regions, sig¬ 
nificant mode-mode interactions should arise: modes that 
extend over large radii should interact with faster growing 
modes localised to small radii. Fluctuations at outer radii 
may be ‘slaved’ to what is going on at smaller radii. Of 
course, distant regions must decouple on some scale, but 
what sets that scale? Also how might all this be connected to 
global field generation and dynamo action? Such questions 


can only be answered by careful numerical simulations. But 
an understanding of the normal mode structure may provide 
useful clues. 

Generalising the equilibrium magnetic field is another 
avenue to explore. Neglected in this paper, azimuthal fields 
are probably the dominant component in a differentially ro¬ 
tating flow; their importance in the disk’s stability could be 
reassessed in an analogous way to here. In addition, a radi¬ 
ally varying magnetic equilibrium should also be revisited. 
Crucially, a vertical field that decays with radius will alter 
the locations of the modes’ turning points; modes of given 
k z and n will extend further outward, changing the nature 
of the disk’s linear response. 

A fourth issue regards the role of large-scale non- 
axisymmetric MRI modes, something we have not touched 
on. Such modes must contend with magnetic resonances and 
consequently, their localisation is more complicated than for 
the axisymmetric MRI (e.g. Curry V Pudritz 1996, Fu & 
Dong 2009). How does this influence, at all, their participa¬ 
tion in disk turbulence? 

Other topics of interest include the breakdown of the 
linear modes due to non-axisymmetric parasitic instabili¬ 
ties (Goodman V Xu 1994), which may provide an alterna¬ 
tive pathway to saturation in some circumstances. Three- 
dimensional cylindrical simulations could probe their be¬ 
haviour and assess their relative importance. A number of 
weakly non-linear analyses of the MRI have been conducted, 
all using local approximations (Umurhan et al. 2007, Jam- 
roz et al. 2008, Vasil 2015). It would be interesting to test 
if analogous calculations are possible in a cylindrical model, 
using the formalism presented in this paper as the (linear) 
starting point. Finally, our cylindrical results could be ex¬ 
tended to vertically structured global disk models, obviously 
making contact with the general theory of Terquem V Pa- 
paloizou (1996) and Ogilvie (1998), but also exploring insta¬ 
bilities in the strong magnetic field limit (cf. Curry V Pu¬ 
dritz 1995, Pessah V Psaltis 2005). The latter may be of par¬ 
ticular interest to the magnetically arrested accretion flows 
around black holes recently simulated (e.g. Tchekhovskoy et 
al. 2011, McKinney et al. 2012). 
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APPENDIX A: MATHEMATICAL 
DERIVATIONS 


A1 Eigenvalue ordering 

In this short section we sketch out a proof showing that 
the eigenvalues e associated with (15) are always less than 
1. For simplicity the boundary conditions are taken to be 
either u r — 0 or d r u r — 0 at r — n, 7*2. The proof for free 
boundaries is a little more involved and we omit its details. 

First multiply (15) by U* and integrate over the do¬ 
main. After integrating by parts and applying the boundary 
conditions the equation can be reworked into 

— 2 _ Jl^l 2 ^ f|drt/| 2 + fr 2 |£/| 2 dr 

f i — 2q \U\ 2 dr + fcf f i — 2q \U\ 2 dr ' 1 1 

We see straightaway that e must be positive. But note also 
that \U\ 2 > r~ 2q \U\ 2 over the entire integration range and 
so the first term in Eq. (Al) must be greater than 1. As a 
consequence, e < 1. 


A2 Eigenvalues in the large k z limit when q — 1 

Here we obtain approximate solutions to the eigenvalue 
equation K u {k z ) = 0 in Section 2.4.2. In the limit of large 
k z both the order and argument of the Bessel function go 
to infinity. We call on the asymptotic expression given in 
Ferreira & Sesma (2008) for the roots of K u (x) when the 
order v is large and imaginary: 

x n « e~™ /2 (v - 2 - 1 / 3 a„e- 27 ri/ V /3 ). (A2) 

Here a n is the n’th root of the Airy function Ai(x). Substi¬ 
tution of v — ik z /e obtains the cubic equation 

e — a n 2 _1/3 /cJ 2/3 e 2/3 — 1 = 0, (A3) 

the correct root of which can be approximated explicitly by 
expanding e in powers of small k^ 2 ^ 3 around 1. 



